function [nd,nf,A,B,reord,CF]=redkw(A,B,CF,nx,lpd,track,tol);
%'hi'

% This function file is the dynamic system reduction program
% as described in "System Reduction and Model Solution Algorithms
% for Singular Linear Rational Expectations Models." 
%
% This version is designed to maximize computation speed in 
% large systems by eliminating operations with large, sparse matrices.
%
% The inputs are the dynamic system elements, i.e., the matrices
% A,B, and CF=[C0 ... Cn], the parameter nx which indicates the number
% of exogenous variables, and the vector lpd, which gives the 
% location of predetermined variables.
%
% Other potential elements are a tolerance levels:
% 'tolp' for determining when elements are nonzero,
% 'tolz' for determining when elements are zero.
%
% The outputs a transformed dynamic system: matrices 
% A,B, PSIF = [PSI0 ... PSIq] and the parameters nf and nd,
% which indicates the number of variables eliminated in the 
% reduction (nf) and the size of the dynamic subsystem
% contained in the specified model (nd).
%
% More formally, the program reduces a singular dynamic
% system of the form:
%
%   A Ey(t+1|t) = B y(t) + C(F) Ex(t|t),
%
% with ny being the dimension of the vector of endogenous
% variables y(t) and nx the dimension of x(t). The notation
% Ey(t+k|t) means the rational expectation of y at date t+k
% given information at t. The result is a dynamic system 
% of the form:
%
%  | 0  0  | |Ef(t+1|t)|   | I     K | |f(t)|    |PSIf(F)|
%  |       | |         | = |         | |    |  + |       |  Ex(t|t)    
%  | 0  I  | |Ed(t+1|t)|   | 0     W | |d(t)|    |PSId(F)| 
%
% with nd the dimension of d(t) and with a reordering of variables
% permitted. Theoretically, the reordering is a (permutation) matrix L, 
% such that  
%               | f(t) |
%  L * y(t) =   |      |
%               | d(t) |
% 
% but this is replaced by an ordering vector in the interest
% of computational efficiency.
%
% The notation reflects the idea that f(t) are variables
% that are similar to "flows" that do not figure in the 
% dynamic subsystem, once they have been "solved out." 
% The remaining variables are intrinsically dynamic and
% hence are called d(t).

% MATLAB versions 4 and 5 differ in some coding.  The
% default is MATLAB 4, but the programs are designed to
% run in either version.  The key differnces between
% the versions in this program are:
%   (a) logical numbers as outputs and indices
%   (b) the length versus max(size(X)) command
 
% set default tolerance for various computations

if (nargin<6)
   track=1;  % tracking is default
end

if (nargin>=7)
   tolz=tol(1);
   tolp=tol(2);
elseif (nargin<7)
   tolz=10^(-6);
   tolp=10^(-10);
   if (track==1)
     disp(' ')
     disp('tolz parameter set to default value of 10^(-6)')
     disp('tolp parameter set to default value of 10^(-8)')
     disp(' ')
   end
end

% determine the size of the y(t) vector

ny=cols(A);

% REORDERING SYSTEM TO PLACE PREDETERMINED VARIABLES LAST

% The KW solution algorithm allows for the variables
% to be reordered, but allows no other transformations
% of variables. The first step is to undertake the intial
% reordering to place the predetermined variables are last.

% The code is designed to run even if there are no 
% predetermined variables (as in Cagan money demand
% model or asset pricing models).

v=version;

if (eval(v(1))>=5)
   npd=length(lpd);
else
   npd=max(size(lpd));   % npd is number of predetermined variables
end

ord=1:1:ny;

if npd> 0   % construct a vector of original ordering
 
      
   % location matrices for predetermined and nonpredetermined variables
   
   npdloc=ones(1,ny);  % initialize location of nonpredetermined variables
   npdloc(lpd)=0;
   pdloc=ones(1,ny)-npdloc;
   
   % construct a reordering vector
   
 if (eval(v(1))>=5)
    reord=[ord(logical(npdloc)) ord(logical(pdloc))];
 else    
    reord=[ord(npdloc) ord(pdloc)]; 
 end
 
% Reorder columns of A and B to reflect reordering of variables.

A=A(:,reord);
B=B(:,reord);

else % If there are are no predetermined variables   
   reord=ord;
end

% LOCATION OF OR INTRODUCTION OF ROWS OF ZEROS IN A

% We begin by specifying an initial number of flows (nf=0).
% As we progressively reduce the system, we always have
% systems of the form

%  | 0  0  | |Ef(t+1|t)|   | I     Nu | |f(t)|    |Cf~(F)|
%  |       | |         | = |          | |    |  + |      |  Ex(t|t)    
%  | 0  a  | |Ed(t+1|t)|   | 0     Be | |d(t)|    |Cd~(F)| 

% The objective is to find more rows of zeros in A and then
% to place more variables into the "f" vector.

ny=rows(A);
nf=0; % number of flows initialized at zero.
nlead=cols(CF)/nx;  
nlead=nlead-1; % take account of C0

% We describe two alternative approaches to introducing rows of zeros 
% The system is partitioned:

% A = | 0  0  |     B = | I     Nu|   C(F) =  |Cf~(F)|     
%     | 0  a  |         | 0      b|           |Cd~(F)|


a=A(nf+1:ny,nf+1:ny);
b=B(nf+1:ny,nf+1:ny);

% We want to get some more zero rows in a, making it

% a  = |  0    0  |
%      | a2lam a2k|


% The standard way of doing this is to compute the singular value
% decomposition of a.  But this may not be necessary if we are
% on the first iteration and there are lots of rows of zeros in A.
% Further, there are usually such rows of zeros because
% most models have some equations with no leads.

% Method 1:

% Check for rows of zeros and then reorder A,B,CF
% based on rows of zeros in A

[sa1 sa2]=sort(max(abs(A')));

% sa1 contains the absolute value of the biggest
% element of each row A.

% sa2 contains indices that sort these from lowest to highest.

A=A(sa2,:);    % rows of A are reordered to place zeros first.
B=B(sa2,:);    % rows of B are similarly reordered.
CF=CF(sa2,:);  % rows of C(F) are similarly reordered.

% determine number of rows of A with effectively zero elements

nz=sum(sa1<tolz);  
ra=ny-nz;
nd=ny;

if (track>0)
   disp(['the initial step of redkw finds ', num2str(nz), ' zero rows in A'])
   disp(' ')
end

if (nz==0)
   
% Method 2: Computation of Singular Value Decomposition

%   Find rows of A that are zero, by computing the svd

[U S V]=svd(a);

% U and V are unitary matrices (U*U'=I and V*V'=I) and S is
% S = |sigma  0 |
%     |  0    0 |
% with sigma being a diagonal matrix with rank(a) postive numbers
% (the singular values) on its diagonal.

% The rank of a matrix is theoretically the number of non-zero
% singular values.  We approximate this here as in the PC-MATLAB
% rank command, where the computation is based on the number of
% singular values greater than a given tolerance level.

s=diag(S);
tol = max(size(a)) * s(1) * eps;
ra = sum(s > tol);

% To utilize the svd, we need to make two system translations

% System Translation #1:

% We now want to transform the system, leaving the flow equations unaltered
% but changing the equations describing the d(t) so that these are 
% in form S*V'*Ed(t+1|t) = U'*b*d(t) + U'*Cd~(F)*Ex(t|t).

UP=U';

A(nf+1:ny,nf+1:ny)=UP*A(nf+1:ny,nf+1:ny);
B(nf+1:ny,nf+1:ny)=UP*B(nf+1:ny,nf+1:ny);

% Translate the C(F) polynomial
CF(nf+1:ny,:)=mconv(UP,CF(nf+1:ny,:),0,nlead);

% System Translation #2:

% We also want to rearrange the ordering of the equations of the system
% so that those with zeros appear before the non-degenerate ones. Formally,
% this involves multiplying by a matrix:
%   | I  0  0|
%   | 0  0  I|
%   | 0  I  0|
% with suitably chosen dimensions of the identity matrices.  This trsansformation
% is necessitated because the singular values are decreasing on
% the diagonal of 'S'.  However, here is accomplished with a row reordering.

rowr=[(1:nf) (nf+ra+1:ny) (nf+1:nf+ra)];

% Translate the matrices
A=A(rowr,:);
B=B(rowr,:);
CF=CF(rowr,:);

if (track>0)
   disp(['the initial step of redkw finds ', num2str(nf), ' A has singular values'])
   disp(' ')
end

end %  singular value decomposition loop

% SYSTEM REDUCTION STEPS

% Intialize iteration 
it=0;

nnf=1; % pretend flows were eliminated on prior iteration

while ((ra<nd)&(nnf>0))
   
% The "d(t) equations" of the system are now in the form:

%   |  0    0  | |Elam(t+1|t)|   | b1lam b1k | |lam(t)|   | PSI1(F)|
%   |          | |           | = |           | |      | + |        | Ex(t|t) 
%   | a2lam a2k| |Ek(t+1|t)  |   | b2lam b2k | |k(t)  |   | PSI2(F)|

% where lam(t) are the non-predetermined d(t) and k(t) are predetermined.
% Notice that this suggests that there are "candidate flows",
% i.e., elements of d(t) attached to behavioral equations that contain
% no leads.  We want to solve for as many of these as possible.

% Get the first nd-max(ra,npd) rows and nd-np columns of the Be matrix, Be11
% rows reorganized.  We take max(ra,npd) since we don't want
% an inconsistent equation system.

ncf=nd-max([ra npd]);
b1lam=B(nf+1:nf+ncf,nf+1:ny-npd);   

% The QR factorization is used here. Matrix b is what we
% are going to decompose. (This b is not the same as 'b' used earlier,
% but rather is a subcomponent of it as explained later.)
% We compute the "QR" factorization of this matrix, i.e., we find matrices
% Q,R, and P such that b P = QR.  Q is unitary, P is a permutation
% matrix and R is upper triangular.  (If the rank of row b is less than full
% then the last rows of R are zero).  

% Hence, if we want to solve the system of equations b x = k, we can do so
% as follows: k = bx <=> k = b P (P'x) <=> Q'k = R (P'x).  Thus, if we are
% willing to reorder the x's (using X=P'x), we can solve for a subset--call
% it x1--as a function of the others: R11 x1 + R12 x2 = k1 => 
% x1 = inv(R11)*{k1-R12*x2}.  In this expression, we have used only the parts
% of the system for which the transformed equations are non-zero.  As a
% consistency check, we also need to verify that the elements k2 =0.

% We are now looking at a system of equations of the form
% 0=b1 * d(t) + c1(F)Ex(t|t), where the 1 denotes that this
% is the subsystem associated with candidate flows, i.e.
% variables to be moved out of the nonpredetermined part
% of d(t) into f(t).  Partitioning the above equation
% further, we have that 0=b1lam*lam(t)+b1k*k(t)+c1(F)Ex(t|t)
% If there are ncf equations, frequently only nnf of 
% these are linearly independent, so we can only solve
% for nnf variables.  The consistency check discussed
% above involves making sure that the "remaining" equations
% (after we have used up nnf components of b1lam) do not
% involve equations such as 0*lam(t)+b1k*k(t)+c1(F)Ex(t|t)
% This involves two ideas.  We must decide whether the 
% coefficient on lam(t) is really zero. If it is and the
% consistency condition is violated, we report whether
% it is violated for k(t) or x(t).

[Q,R,E]=qr(b1lam,0);
% System Translation #3:

reords=reord(nf+1:ny-npd);  %locations of lam(t)

reord=[reord(1:nf) reords(E) reord(ny-npd+1:ny)];

reordl=[1:nf E+nf (ny-npd+1):ny];

A=A(:,reordl);
B=B(:,reordl);

% This reordering implies that the ncf=nd-max(ra,npd) candidate flows
% are ordered first in the d(t).

% System Translation #4:
% multiply the ncf equations of the "candidate flows" by Q'

B(nf+1:nf+ncf,nf+1:ny)=Q'*B(nf+1:nf+ncf,nf+1:ny);

CF(nf+1:nf+ncf,:)=mconv(Q',CF(nf+1:nf+ncf,:),0,nlead);

% Solving for the candidate flows.

% We determine how many linearly independent rows of the equations there are
% by finding the total number of elements in each row that are close to zero;
% computing the number of rows with all zeros and subtracting from the number
% of rows.

r=abs(diag(R(:,1:rows(R))));  
tol = max(size(b1lam)) * r(1) * eps;  
nnf = sum(r > tol);

if (nnf>0)
% We check the elements of the linear system if nnf < ncf
if (nnf<ncf)
   qx=max(max(abs(CF(nf+nnf+1:nf+ncf,:))));
   if (qx>tolp)
     clc
     disp('The redkw program--system reduction--has encountered')
     disp('major problems with flow equations and exogenous variables')
          disp('This can arise for either of two reasons. First, you')
     disp('may have a well-specified model containing an equation')
     disp('of the form:  0=sm*lamj(t)+c1(F)*Ex(t|t),')
     disp('where lamj(t) is nonpredetermined and sm is sufficiently small')
     disp('that elements of c1(F)/sm falls above the "tolp" value, which is')
     disp(tolp)
     disp('To explore this run the debugkw.m program')
     disp('  ')
     disp('strike any key to continue')
     pause
     clc
     disp('The other possibility is that your model is not')
     disp('well-specified, since KW prove that a model for which')
     disp('any solution exists cannot fail in this way.')
     disp('To see an example of how this can fail, look at')
     disp('crash examples on the replication diskette.')
     disp(' ')
     disp('strike any key to continue')     
    pause
  end
  qp=max(max(abs(B(ny-npd+1:ny,:))));
  if (qp>tolp)
     clc
     disp('The redkw program--system reduction--has encountered')
     disp('major problems with flow equations and predetermined variables')
     disp('This can arise for either of two reasons. First, you')
     disp('may have a well-specified model containing an equation')
     disp('of the form:  0=sm*lamj(t)+bkj*k(t),')   
     disp('where lamj(t) is nonpredetermined and sm is sufficiently small')
     disp('that elements of bkj/sm falls above the "tolp" value, which is')
     disp(tolp)
     disp('To explore this run the debugkw.m program')
     disp('  ')
     disp('strike any key to continue')
     clc
     disp('  ')
     disp('The other possibility is that your model is not')
     disp('well-specified, since KW prove that a model for which')
     disp('any solution exists cannot fail in this way.')
     disp('To see an example of how this can fail, look at')
     disp('crash examples on the replication diskette.')
     disp(' ')
     disp('strike any key to continue')     
    pause
  end
end

% new number of flows: nf+nnf

nf=nf+nnf;
nd=ny-nf;

% undertake classical system reduction

[A,B,CF]=csr(A,B,CF,nf,nnf,nlead);

nlead=nlead+1;

it=it+1;

if (track>0)
   disp('   ')
   disp(['after iteration # ', num2str(it), ' of redkw.m'])
   disp(['number of orginal variables (ny) = ', num2str(ny)])
   disp(['number of flow    variables (nf) = ', num2str(nf)])
   disp(['number of dynamic variables (nd) = ', num2str(nd)])
end

%  USE THE SVD TO COMPUTE RANK AND PREPARE FOR NEXT ITERATION

a=A(nf+1:ny,nf+1:ny);
b=B(nf+1:ny,nf+1:ny);

[U S V]=svd(a);
s=diag(S);
tol = max(size(a)) * s(1) * eps;
ra = sum(s > tol);

if (ra<ny-nf)  % additional reduction steps must be undertaken
   
UP=U';

A(nf+1:ny,nf+1:ny)=UP*A(nf+1:ny,nf+1:ny);
B(nf+1:ny,nf+1:ny)=UP*B(nf+1:ny,nf+1:ny);

CF(nf+1:ny,:)=mconv(UP,CF(nf+1:ny,:),0,nlead);

rowr=[(1:nf) (nf+ra+1:ny) (nf+1:nf+ra)];

A=A(rowr,:);
B=B(rowr,:);
CF=CF(rowr,:);

end %  singular value decomposition loop

else % nnf=0
   disp('  ')
   disp('There are no new flows that can be isolated')
   disp('(nnf=0) and the matrix A is singular.')
   disp('The system cannot be reduced.')
   disp('  ')
   disp('run the debugkw.m program')
   disp('also see crash examples')
   disp('  ')
   disp('The matrix b1lam is:')
   disp(b1lam)
   disp('  ')
   disp('strike any key to continue')
   pause
   disp('  ')   
end

end

% OUTPUT OF REDUCED SYSTEM

if (ra<ny-nf)
  disp('dynamic system contains irreducible singularity in A')
else
   
   % Put the system in "standard form", by inverting "a"

   %  | 0  0  | |Ef(t+1|t)|   | I     K | |f(t)|    |PSIf~(F)|
   %  |       | |         | = |          ||    | +  |        |  Ex(t|t)    
   %  | 0  I  | |Ed(t+1|t)|   | 0     W  ||d(t)|    |PSId~(F)| 
   
   
  B(nf+1:ny,nf+1:ny)=A(nf+1:ny,nf+1:ny)\B(nf+1:ny,nf+1:ny);
  ia=inv(A(nf+1:ny,nf+1:ny));
  A(nf+1:ny,nf+1:ny)=ia*A(nf+1:ny,nf+1:ny);
  CF(nf+1:ny,:)=mconv(ia,CF(nf+1:ny,:),0,nlead);
end

if (track>0)
   disp(' ')
   disp('system reduction completed')
   disp(' ')
end

nd=ny-nf;



